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Abstract 

To investigate order-order interfaces, we perform multimagnetical 
(U i Monte Carlo simulations of the 2D and 3-D Ising model. Following 



"^ ' Binder we extract the interfacial free energy from the infinite volume 

limit of the magnetic probability density. Stringent tests of the nu- 

^ , merical methods are performed by reproducing with high precision 

^ [ exact 2D results. In the physically more interesting 3D case we esti- 

mate the amplitude Fq of the critical interfacial tension F^ = F^f^ to 
be Fq = 1.52 zb 0.05. This result is in good agreement with a previ- 
ous MC calculation by Mon, as well as with experimental results for 
related amplitude ratios. In addition, we study in some details the 
shape of the magnetic probability density for temperatures below the 
Curie point. 
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1 Introduction 

Over the last fifty years, and presumably longer, there has been continuous 
interest in the properties of interfaces in 2D and 3D Ising models. In the 2D 
case a number of remarkable exact results could be derived |I|, |], ^ , whereas 
in 3D Monte Carlo (MC) simulations 0, |^, |^, |^ have become the dominant 
method. Of particular interest are amplitude ratios p, ^ which involve 
the amplitude Fq of the critical 3D interfacial tension, because they can be 



compared with a variety of experimental results [|r^, ^ for fluids which are 
supposed to populate the Ising model universality class. Additional interest 
in Fq is created by the fact that it controls the non-universality of the critical 
wetting transition, a subject in which one has seen continuous interest over 
the last years |1T2 |. 



Let us introduce the notation. Spins Sj = ±1 are defined on sites of a 
square lattice of volume V = L^ with periodic boundary conditions and the 
symbol < i,j > is used to denote nearest neighbours. The partition function 
of the D-dimensional Ising model is given by 

Z = Y. exp{-pH). (1) 

configurations 

The Ising Hamiltonian H 

H = Hi - hM with Hj = - J2 ^i^j ^^^ ^ = J2 ^i (2) 

contains the nearest neighbour interaction term Hj and a term which couples 
the external magnetic field h to the magnetization. We denote the magnetic 
probability density by Pl{M). Often it is preferable to use the magnetization 
density m = M/L^ and 

PLim) = Pl{M). (3) 

When convenient arguments are suppressed altogether. Let us consider the 
case of a vanishing magnetic field /i = 0. In the broken region j3 > (3c 
the magnetic probability densities Pl{M) are double peaked. We denote the 
positions of the maxima by ±M™^^. As the model is globally Z[2) symmetric, 
we exploit the symmetry PiiM) = Pl{-M), define P^'^^ = Pl(M]^^^) = 
Pl(-M™''''), and confine our discussion to M > 0. For Ml < M"^""" the 
probability density takes then its minimum at P™™ = Pi(0). 



The surface tension F'^ is defined as the the free energy per unit area of the 
interface between coexisting phases. Here we consider order-order interfaces 
in the broken region j3 > Pc- The asymptotic behaviour for /3 — ;> /3c is 

F' = F^t^" + F^tf'+\ (4) 

where t = (1 — Pc/P) is the reduced temperature. Widom's scahng relation 

m 

/i = (D-l)z/, (5) 

connects the exponent /i with the critical exponent z/ of the correlation length. 
In 2D n = u = 1, together with Fq = 1.762... and F^ = 1.216... follows from 
exact results |jl|. In 3D approximate methods have to be used, and z/ = 0.630 



is in agreement with [p!^ , |I^ as well as consistent with many more estimates 
of the extensive literature. Therefore, we use /i = 1.26 in 3D throughout this 
paper. When combined with other amplitudes, F^ enters into a number of 
universal ratios |^, ^ . The other amplitudes are fairly accurately known, see 



for instance ||T6[, and the uncertainty in the ratios stems mainly from F^. 

Numerical calculations of interface tensions have remained subtle and a 
large number of methods have been tried. Old estimates of Fq |Q led to 
amplitude ratios which are in disagreement with the experimental results for 
fluids Ijl^, 0, whereas Mon obtained consistency by calculating the ex- 



cess energy between a system with an interface imposed by an antiperiodic 
boundary condition and a system without an interface with periodic bound- 
ary conditions. Another promising approach is to calculate the interfacial 
tension from correlation functions |l^, 0. In this paper we turn back to 
the approach of Binder |0, ^ and combine it with a new simulation method 

mm- 

The interface tension F^ can be defined from the finite size scaling (FSS) 
behaviour of pmm/pmax^ which for large L takes the form 



P, 



min 
L 



pniax 



const LP exp -2L^-^F^(1 + aiL'^ + aaL"' + 0{L-^)) (6) 



on lattices with periodic boundary conditions. Physically this definition as- 
sumes that at values of the mean magnetization M = a rectangular domain, 
enclosed by two interfaces and spanning the lattice via the periodic bound- 
ary condition, is formed. Interface shapes other than planar, see for instance 



Figure 7 in section 6, are due to capillary waves which have to be included 
in the definition of F^. It is convenient to introduce the function 

pniax 

^z.(M) = L-(^~i) IniTTiTV ^nd M^) = ^l{M). (7) 

The finite lattice interfacial tension is then defined as 

F[ = 1^l{0) = ^MO), (8) 

and F^ = linii^oo F[. It follows that the infinite volume limit 

ilj{m) = lim ipL{m) (9) 

exists and takes values in the range [0, F^]. 

Originally, the calculation of F^ via (^ suffered from a problem of prin- 
ciple. To get a representative sample of interfaces one has to generate many 
configurations with M = on large lattices. But, in a canonical simulation 
this is precisely impossible because of (H). This equation means that the rel- 
evant configurations are exponentially suppressed, implying a slowing down 
of the simulation ~ exp(-|-2L^~^F'*). Recently this difficulty was overcome 
by the proposal to perform the MC simulation for a multimagnetical ensem- 
ble [0, a natural extension of the multicanonical ensemble introduced in 
1^ . The slowing down becomes then reduced to a power law close to ~ \^^. 



Details of this approach are given in the next section. 

In section 3 the numerical methods are introduced which we employ to 
cope with equation (^ on the basis of our multimagnetical MC data. Section 
4 summarizes our 2D and section 5 our SD results. Besides the surface ten- 
sion, this includes tunneling times and data for the magnetic susceptibility. 
The 2D Ising model serves essentially as a testing ground for the numerical 
approach, whereas the 3D results are of more physical interest. Section 6 
discusses the shape of the magnetical probability density in greater details. 
A brief summary is given in our final section 7. 



2 The multimagnetical ensemble 

The standard MC algorithm would only sample configurations correspond- 
ing to P™° if one could generate of the order 1/P™° or more statistically 



independent configurations. We overcome this difficulty by sampling config- 
urations with a multimagnetical weight factor 

P™'"(M) ~ n{M)-^ for < M < M]^^^. (10) 

Here n{M) is the magnetical density of of states (fixed temperature (3~^). 
For any reasonable approximation the resulting multimagnetic probability 
density will be approximately flat. For reasons of physical intuition the 
parameterization 

p™'»(M) = exp(ai(M) + /ii(M)/3M - (3Hi) for 0<M<V (11) 

is useful. For constant a, h the weight factor exp(a + hf3M — j3H) describes 
a canonical ensemble in an external magnetic fleld h. The function /^^(M) 
deflnes an M-dependent effective external magnetic field such that the re- 
sulting multimagnetical probability density is fiat for M < M™^^. As M 
changes in steps of two, the definition of h^lM) for 2 < M < Mf^^^ — 2 is 

exp [2/3/..(M)] = ^^^^. (12) 

Including now M = and M > M™'^'' leads to 

JO for M = and M > M"'^^; 
hL{M) - [ 2~^f]-^L^-^[^i,^(M)-^L{M + 2)] elsewhere {M>0). ^^' 

For M < 0, hilM) is defined through hi{M) = —hi{—M). The purpose 
of the function aL{M) is to ensure continuity of the multimagnetical weight 
factor. This leads to the recursion relation 

aL{M + 2) = aL{M)+P[hL{M)-hL{M + 2)] (M + 2), ^^(0) = 0. (14) 

Once hiiM) is given, aiiM) follows automatically. The standard Markov 
process is well-suited to generate configurations which are in equilibrium with 
respect to this multimagnetical distribution. The h = canonical probability 
density Pl{M) is obtained from PJ^'^{M) by: 

Pl{M) = cPr'^iM) exp{-aL{M)-hL{M)f3M). (15) 

The constant c is obtained by imposing the appropriate normalization on 
Pl{M), which we choose P^^^ = 1 for the purposes of this paper. 
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How do we get now the inultiinagnetical function h{M)l For small sys- 
tems Pl{M) can accurately be calculated by performing standard heatbath 
or Metropolis simulations, and h^lM) follows directly from (0JT^). On larger 
systems we get h^lM) by making every time a FSS prediction of \I/2,(M) from 
the already controlled smaller systems. To optimize the parameters we may 
perform a second run and in some cases further runs to increase our statistics. 
We now explain details of the FSS prediction. The outlined method is simple 
and worked sufficiently well for our practical purposes, but the approach is 
by no means unique and can likely be improved. 

From the definition of the interface tension we know that iPl{0) scales 
like 



aO 



MO) = ^F' + JD^^ (16) 

where we neglect ln(L) terms and a° is an unknown constant. Further, 



I I L T 



j^max^^D jg ^ function of the lattice size and scales like 



^max _ ^max / ^ ^ _^_ \ _ (^7) 

We employ this equation for re-adjusting the magnetization m^ in general, 
and use a scaling law of the form (ITO) for the entire function ?/'j^(m): 



^LirriL) = V^(m) + j^, (18) 

where iplm) and a™ are two unknown, m-dependent constants. To calcu- 
late the multimagnetical parameters for a new lattice size we just take our 
two largest lattices, start from the m^ on the smaller lattice, determine the 
readjusted magnetizations on the larger lattices and compute for this set of 
points ipL{m) for the new lattice. A linear interpolation between these points 
gives us our final approximation of ipLlm). Using \E'(M) = ipim), the new 
multimagnetic function hi{m) follows again from ([TB|). It may be noted that 
due to the choice of our normalization a™ = at rn^^^ for all L. Since 
we only need an approximate estimate of hiiM), it seems to be harmless 
that \E'l(M) is a noisy estimate. It is our experience that the noise averages 
out over sufficient distances in m. In Figure 1 we give an example of an 
un-reweighted distribution which we get by this method. 



As before |18|, |I9| we define the tunneling time tl as the average number 



of updates needed to get from a configuration with M = M™™ = —M 
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max 



to a configuration witli M = Mf^^^ and back. Due to (|T0|) we expect that 
the standard Markov process leads to a one dimensional random walk like 
behaviour of the magnetization M. Because of M™^^ ~ V the slowing down 
becomes then reduced to a V"^ power law. Our data in the subsequent sec- 
tions confirm this up to a small correction, which indicates that the random 
walk behaviour is not entirely perfect. The deeper we are in the broken re- 
gion, the more more dramatic is the improvement which we expect. Because 
F"* increases with /?, this follows from (H). 



3 Numerical analysis of the surface tension 

To calculate F£ via (^§) requires the maxima and the minimum of the 
probability density Pl{M). To determine the minimum causes no problem 
because we know p™™ = P^(0). To estimate P™^^ we exploit the exactly 
known symmetry Pl{M) = Pl{—M): P™^^ is estimated by looking for the 
absolute maximum of the distribution on one side, taking the reflected value 
on the other side, and averaging over both possibilities. 

The final estimate of F'^ requires a FSS extrapolation towards L = cxd. 
Equation (|^) converts into the fit 

Here c/L^ takes care of the L~^ corrections to F^ in (^. However, the 
accuracy of our F[ data will not allow a four parameter fit. The obvious three 
parameter fit (HIi) is obtained by setting c = 0. Further, it is numerically 
difficult to include the a/L^~^ and the bln{L)/L term. Therefore, the three 
parameter fits defined by 6 = (III2) and a = (III3) are also considered. 
Corresponding two parameter fits (II2) and (II3) follow by demanding c = 
in addition. Although interactions (correlations) between the two interfaces 
fall off exponentially with their distance, and hence the system size, for small 
systems they may be more important then the L~^ correction to F'^ in (|^). 
Consequently one may try the fit 

Fl = F. + ^ + !!J^ + ce-- (20) 



with either a = or 6 = 0. 

One hkes to include data from as many lattices as possible,while keeping 
the fit self-consistent. The systematic approach towards this fitting problem 
begins with fitting the data from all the lattices and monitoring x^ of this fit. 
The fit is accepted if the obtained x^ value results into a reasonable likelihood 
(here we choose > 10 % ) that the discrepancies between the fit and the data 
points are due to statistical fluctuations. If this is not the case, the smallest 
lattice is supposed to suffer from strong finite size corrections of a form not 
included in the fit. Omitting the smallest lattice, the fit is carried out again. 
The whole procedure is repeated until an acceptable, called self-consistent, 
fit is obtained (or until one runs out of sufficiently many data points). 

The statistical error of F£ has a somewhat peculiar finite size behaviour. 
According to our discussion of the previous section we assume that the multi- 
canonical simulation slows down like L^'^'^'^^^ with i] > small. In this paper 
we give the statistics in sweeps. One sweeps updates each spin on the lattice 
once. For a constant number of MC sweeps we have 

A pmin 

t±Lk L^l+^)D/2 (21) 

pinin ' \ ' 

because the number of independent entries of -Pl(O) decreases ~ L^^^+^^^ . 
Now (remember our normalization p™'^^ = 1) 



-In 

and, consequently. 






L 



pmiri 



A pmin 

2i:P-i)p« ± ^±L^ (22) 

pmin V / 



[/\Flf ~ L^n-i)D+2 ^ ^2-D_ (23) 



This equation implies that in D = 2 a constant number of sweeps keeps the 
Fl error bar approximately constant, and in D = 3 this may even be achieved 
with a decreasing statistics (~ 1/a/L for r^ = 0). Of course, the statistics has 
to be chosen such that on the largest lattice still a sufficiently large number of 
tunneling events {i.e independent data) is obtained. Otherwise a meaningful 
statistical analysis is impossible. Twenty or more tunneling events is a a 
good number for practical purposes and numbers as small as four do still 
allow some kind of reasonable analysis |pO[| . 



We pursue two ways to estimate the amplitude Fq defined by (H). For a 
sufficiently small reduced temperature 

Fo^i := FH-^" (24) 

approximates Fq well. On the other hand, when we keep the data time series 
re- weighting techniques [^, ^ allow to calculate -F**(/?) in a neighbourhood 



of the simulation point. This allows another estimate of F^: 

F,, :^ M,Pt With n^§^ ^ '%:i:T - <^^' 

Estimates (p^) and (^) differ by the next term in the expansion (1): 

Fl = f,r\F^,-F^,) (26) 

In principle one may use this to estimate Ff and, subsequently, use (|^) to 



improve the estimate (plf) of F^. In practice the estimate F^2 provides only 
a consistency check, because Fq2 turns out to be too noisy to compete with 



-^01- 



4 2D estimates 

For the two dimensional Ising model we performed multimagnetical simu- 
lations at the critical temperature (3 = (3c = ln(l + v^)/2 = 0.44068..., at 
(3 = 0.47 and 0.5. The results are collected in Tables 1-3. In each run 
additional 200, 000 initial sweeps without measurements were performed for 
reaching equilibrium with respect to the multimagnetical distribution. To 
give an example. Figure 2 depicts the (3 = 0.5 magnetic probability densities 
Piijn) for the range L = 10 — 100. Beyond F£, F£ and t^, as introduced in 
the previous sections, another quantity included in Tables 1-3 is the magnetic 
susceptibility 

XL = V^w?pL{m). (27) 



In the large L limit its leading term will approach the exactly known p2 
mean magnetization squared: {imPY . Recent progress in FSS theory raises 



interest |^ in the subleading terms, which even in 2D are not known ana- 
lytically. 

For all f3 values self-consistent three parameter fits (JT^IIIi) are possible. 
The obtained parameters b are compatible with zero, and the appropriate 
two parameter fits ([T9I-II2) do not decrease the lattice range over which self- 
consistency is achieved: L = 14 - 100 (/3 = /?J, L = 20 - 100 {(3 = 0.47), 
and L = 30 — 1002 {j3 = 0.5). Including now the L~^ term, i.e. turning to 
fit ([T9I-III2), increases the self-consistency range considerably: L = 4 — 100 
{/3 = (3c), L = 10- 100 {(3 = 0.47), and L = 8 - 100 {(3 = 0.5). Together 
with the F[ data these fits are depicted in Figure 3, and we rely on them for 
our final F^ estimates. Exact results follow from Onsager's [|l| equation 



2/3 -In 



1 + e-^f 
1 _ e-2/3 



(/3 > (3c). (28) 



In Tables 1-3 we have included our L = 00 estimates as well as the exact 
values. Good agreement is found in all cases. 

To calculate the amplitude Fq we employ equations (^4]) and (^51). We 
list in Tables 1-3 the corresponding numerical estimates F^^ and Fq2 together 
with their exact values which follow from Onsager's formula (^). It is seen 
that Fq2 is systematically larger than F^^. According to (^) this implies 
Ff > 1. Equation (^) defines then a third, improved approximation Fq^ 
of Fq which is even smaller than Fq^. This may be illustrated with exact 
results: Fq^ = 1.76 for (3 = 0.47 and Fgg = 1.75 for /3 = 0.5 are obtained. 
Both values are good approximations of the exact results Fq = 1.74. The 
slight improvement for (3 = 0.5 is accidentally caused by next order terms. 
Unfortunately, the numerical estimates of Fq2 from the F[ turn out to be 
too noisy to make this a feasible practical method. Still, the increase of the 
Fq^ estimates with t as well as the Fq2 give correctly the trend. In fact the 
direct two parameter fit (|^) of our two rather accurate F* estimates give 
Fq = 1.765 ± 0.003. In conclusion, the employed methods pass well testing 
versus the exactly known 2D Ising model results. 

To estimate the slowing down of the multimagnetical simulations we fit 
the Ti values of Tables 1-3 to the form tl ~ L^+^^, corresponding to the 
suspected slowing down ~ V"^^"^ in updates. We find r] = 0.012 {(3 = (3c) , 
rj = 0.045 (/? = 0.47), and r] = 0.048 (/3 = 0.5). The actual improvement is, 
of course, most impressive for /3 = 0.5, as the exponential slowing down of 



10 



the canonical simulation is there strongest. See [|T^] for a detailed comparison 
with standard heat bath simulations at this f3 value. At /3 = /3c the standard 
simulation does not slow down anymore exponentially. Still we find improve- 
ment by a constant factor ^ 3, and this situation is depicted in Figure 4. 



5 3D estimates 

We performed simulations at /5 = 0.227, 0.232 and 0.2439. The corresponding 
temperatures are all well above the roughening temperature /^^ ~ 1.85 /3c ~ 
0.41 |2|]. Our results are collected in Tables 4-6. As in the 2D case 200,000 



additional, initial sweeps were performed in each run for reaching equilibrium. 
Again, we give an example for the magnetic probability densities. Figure 5 
illustrates (3 = 0.2439 for L = 8 - 32. 

In Figure 6 we display the effective tensions as functions of the lattice size 
together with asymptotic fits. We notice that finite size effects play a more 
important role in three than two dimensions and are more complicated, too. 
The non-monotone behaviour shows that it is necessary to use large enough 
lattices to estimate the interface tension. For each (3 value we could only 
generate few lattices of sufficiently large size. Consequently, we find that our 
F[ data are only well suited for two parameter fits. In contrast to 2D, it turns 
out that the ln(L)/L term seems to be relevant. The two parameter fit ([T9| - 
II3) is preferred by the data when compared with (^^112), presumably both 
terms would finally contribute. The self-consistency ranges are L = 24 — 32 
{(3 = 0.227), L = 20 - 32 (/3 = 0.232), and L = 24 - 32 (/3 = 0.2439). We 
include the thus obtained estimates in Tables 4-6. The accuracy of our data 
does not allow to determine the exponent p in (^ reliably. Obviously, the 
situation is far less satisfactory than in 2D. Larger lattices are desirable, 
but could not be simulated with our present computational resources. To 
get a handle on the order of magnitude of the expected systematic errors, we 
tried to extend the self-consistency ranges by including further fit parameters 
(|19D, (^). With one exception these fits are either unstable or have a bad 
goodness of fit. The exception is the (3 = 0.232 point. The fit ( pO|) with a = 
leads to the self-consistency range L = 6 — 32 and F* = 0.03044 ± 0.00025. 
By comparison with the results reported in Table 5, one may argue that the 
results of Tables 4-6 are presumably afflicted by systematic errors of 2% to 
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4%, thus exceeding the statistical errors. 

The value j3 = 0.232 was chosen, because it enables a comparison. In 
cluster improved estimators were used to calculate correlations for L = 8 — 14 
in a cylindrical geometry. Fitting the obtained tunneling mass gaps yields 
-F* = 0.03034 ± 0.00015 for the surface tension. This is in remarkably good 
agreement with our fit pO| ) of the range L = 6 — 32, although this might be 
accidental. In any case, one seems to be on the save side by allowing for the 
discussed systematic error. All our F^ values are much higher that the old 
estimates of |^, which had to rely on far too small lattices. 

To estimate the amplitude Fq we first calculate Fq^ ( p^ and Fq2 (p5D, as 
in the previous section, and include the results in Tables 4-6. In contrast 
to the 2D case, no clear trend of Fq^ with increasing t is observed. In view 
of this, we tend to attribute their discrepancies to be due to the systematic 
errors in taking the asymptotic limits limi_^oo -^l (/^) • Consequently, a fair 
estimate is obtained by averaging over our three values and allowing for about 
3% systematic errors. This gives 

F^ = 1.52 ±0.05. (29) 

Surprisingly, our 3D results for Fq2 suffer less from statistical noise than 
in 2D. As before the Fq2 ^^^ systematically higher than the corresponding 
Fq]^. This indicates that a further downward correction of (pOf ) might be 
necessary. However, presently involved ambiguities are too severe to dare 
such a correction. Instead one should first try to push the multimagnetical 
3D Ising model calculations towards larger lattices, an enterprise which will 
become feasible with the next generation(s) of fast workstations. 

Again, we perform power law fits of the tl data to estimate the slowing 
down of the multimagnetical simulations. We fit the t^ of Tables 4-6 to the 
form Tl ~ L^~^^^, corresponding to the suspected slowing down ~ V'^'^'^ in up- 
dates. We find r] = 0.001 (/3 = 0.227), r] = 0.006 {(3 = 0.232), and r] = 0.030 
(/9 = 0.2439). The real improvement is most dramatic for (3 = 0.2439. For 
the largest lattice lattice (L = 32) it is about exp(2 x 0.074 x 32^) ^ 10^^ 
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6 The shape of the probabihty distributions 

The behavior of the magnetic probabihty density at smah values of the mag- 
netization is related to the existence and properties of droplets and domains 
of one of the pure phases floating in a medium of the opposite phase. On a 
given lattice properties of domains will strongly depend on the chosen bound- 
ary conditions, which in our case are periodic. To proceed we shortly review 
2D results of 0, ^. The analytical form of i'im), defined by @, can be 
given in terms of three temperature dependent parameters m^, w and m'^: 



w y/m^ — m'^ = const = 2F* for — m'^ < m < Trf] 

HM) = { (30) 

w JmP— I m I for | m |> m^. 

There exists a finite interval in the absolute value of the magnetization, where 
the function ijj is constant, namely twice the value of the interface tension 
along one of the main axes of the lattice. At the value oim = Trf the function 
tp develops a cusp. Following a square root behavior, ip then approaches zero 
at the value of the mean magnetization m^. 

Let us outline the qualitative physical picture behind the above formula. 
At magnetization m k. Q half of the spins belong to a domain with a char- 
acteristic magnetization corresponding to one of the pure phases, e.g. to the 
state characterized by -l-m". Consequently the other half of spins are con- 
tained in a domain structure with typical magnetization —m^. In between 
these domains a surface is formed, and dominant contributions to the parti- 
tion function at m ~ arise from configurations with a minimal excess free 
surface energy. These are configurations where a domain strip is closed via 
the periodic boundary conditions of the lattice. In Figure 7 we display a 
typical configuration contributing at m ^ on the 100^ lattice. One clearly 
observes a strip with sizable fluctuations added. Such configurations allow 
translation of one interface without changing the energy or the number of 
states. As the magnetization is changed under such a translation one finds a 
constant behavior of ■?/' in a finite interval | ?7i |< nf. At m = nf yet another 
mechanism takes place. It becomes more favorable for the system to form a 
droplet - or bubble - of one pure phase floating in a sea of the opposite phase. 
This happens when the surface energy of a droplet is smaller than the surface 
energy of a domain. The system enters a region of configuration space where 
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single or multiple droplets dominate. Figures of such configurations are given 
in 
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These statements can be made quantitative for the 2D Ising model. Let 
us consider the point (3 = 0.5. Onsager's [|l| interface tension determines 
w ^/mF^^mP = 2F' = 0.45612 for (3 = 0.5, and due to Yang [|2] m° = 0.9113 



for /5 = 0.5. The remaining free parameter is the cusp value of the magne- 



tization m^. Its value follows from the Wulf construction |^ of classical 
equilibrium shapes of droplets. Relying on analytic results of Rottmann and 
Wortis we display in Figure 8 equilibrium shapes of droplets for several 
temperatures below Tc. At very low temperature a quadratic shape is as- 
sumed with small rounding effects at its corners (at zero temperature the 
shape is exactly quadratic). Closer to the transition point the droplet shapes 
take a circular form, while exactly at the transition point due to rotational 
invariance the droplet shape is a circle. At jS = 0.5 the droplet shape ex- 
hibits only tiny deviations from a circular shape. When the surface energy of 
a classical droplet equals the surface energy of the rectangular domain 2LF^, 
the transition from the droplet phase to a percolated state takes place. For 
a circular droplet shape one finds the relation 

"" " "" ll-4/(3v/^) for D = 3, ^"^^^ 

while for a quadratic shape m'^ = m°/2, and m'^ = m^/3 for a cubic shape. 

In Figure 2 we already presented our numerical /3 = 0.5 probability distri- 
butions piiiTT-) on lattices ranging from 10^ to 100^. As predicted by theory 
we observe with increasing lattice size the unfolding of a flat region at val- 
ues of the magnetization close to zero. In order to compare the theoretical 
scenario with our simulation, we construct estimators of the function ipi by 
performing three parameter fits of our numerical data with the analytical 
form (^0|). The results for the fitted parameters m'^ are displayed in Figure 
9 as a function of the inverse linear lattice size 1/L. We observe a rather 
rapid variation of the estimate of the cusp magnetization value, indicating 
that even for our largest lattices the estimated function ipi is relatively far 
away from its asymptotic value. Nevertheless, the data do allow a consistent 
extrapolation to the theoretical value as obtained from (^) (dashed curve in 
the plot). 

In Figure 10 we display the extrapolation of our finite volume ipi data 
to their infinite volume limit. Similarly as in section 2 we assumed the 
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leading finite size correction to be proportional to 1/ L. Excluding lattices 
of size smaller than 60 from the fit, we obtain our infinite volume data. It 
is obvious from our discussion in section 2 that a successful extrapolation of 
this type determines the multimagnetical parameters. We have fitted these 
infinite volume data with the analytical form of equation (0), where we left 
the parameters rrf, w free and fixed m° to its exact value. We obtain a good 
fit consistent with the predicted theoretical shape. The fitted value of the 
cusp magnetization comes out to be m'^ = 0.32(1), which is to be compared 
with the theoretically predicted value m"^ = 0.33. 

In Figure 5 we displayed the 3^ probability densities PLijn) at /? = 
0.2439. We also find indications of a flat region. However our lattices are 
too small to allow an estimate of the ?)D cusp magnetization by the above 
method. Therefore, we tried another approach, which is on less solid theo- 
retical grounds, but seems to work also for the 3D Ising model, provided (3 is 
sufficiently large. In this approach we try to determine the cusp magnetiza- 
tion by measuring the total interface area as function of the magnetization. 
This area is just the number of nearest-neighbours on the lattice with oppo- 
site spin. We assume that this area is either dominated by droplet (s) or, after 
percolation, by the rectangular domain. For D = 3, /5 = 0.2439 and various 
L we plot in Figure 11 the surface area, divided by the number of links. 
For large enough L one can easily distinguish three different parts in these 
functions. In the vicinity of of the negative mean magnetization the surface 
area increases rapidly: more and more small droplets are created. Then the 
droplets join to larger bubbles and these bubbles grow in size. The increase 
of the total surface area becomes then somewhat irregular, as it may de- 
crease by joining droplets, but over-all the increasing trend continues. After 
reaching a final maximum value, the total interface area decreases somewhat, 
becomes then fairly fiat (very gently decreasing) until m = is reached. For 
positive magnetization m the picture is just refiected. The final maximum 
is obviously associated with percolation of the largest bubble and the subse- 
quent decrease has to be an entropy effect. We argue that there is a range in 
the magnetization where the rectangular domain is energetically favourable, 
but the free energy is still minimized by configurations containing also some 
droplets. With further decreasing magnetization the energy gap between the 
droplets and the domain grows and finally the domain dominates. 

To estimate the cusp magnetization nf, we take the value of the sur- 
face area at ttt, = (where one rectangular domain dominates most), and 
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determine the point ml of the increasing part of the function with identi- 
cal surface are. This is the point where the droplets and bubbles have the 
same total surface area as the domain. For those systems where this analysis 
is possible the ml values are given in Table 7. We cannot determine m^ 
for all our systems, since the time data series is needed, which we did not 
keep on disk for each case. We notice that the cusp magnetization m.'i as 
function of the lattice size is almost constant. Allowing a 1/L^^~^) correc- 
tion, we extrapolate to L — ;► cxd. We find m'^ = 0.35 ± 0.02 for the 2D Ising 
model at f3 = 0.5. In 3D we find m"^ = 0.33 ± 0.01 at f3 = 0.2439. This is 
in good agreement with the value one obtains using ( pi] ) and our estimate 
m^ = 0.709 of the mean magnetization. As all our /? values are well below 
the roughening point, one expect a circular shape and hence (^Tj) to be valid. 



7 Summary 

Multicanonical simulations allow to study the magnetic probability density 
in the broken region with a hitherto unreached precision. Our numerical 
calculations for the 2D Ising model agree well with exact results. For the 
3D Ising model we obtain new surface tension estimates. At /3 = 0.232 good 
agreement is found with and our amplitude estimate agrees with Mon [^]. 
When presenting the results of this paper at the Jiilich workshop "Dynam- 
ics of First Order Phase Transitions" EHl, we became aware of two other 



recent investigations [27, 28 1 which employ different methods to address sim- 



ilar problems as the present paper. Our results are in good agreement with 



and seem also to be consistent with 28 
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Tables 



Table 1: The tensions -F£ as function of L in D 



2for/5 = /3, (t = 0). 



L 


No of sweeps 


F[ 


Fl 


tl 


XL 


2 


4 


10^ 


0.53490 (42) 




16.602 (1) 


3.4974 (9) 


4 


4 


10^ 


0.41288 (67) 




90.58 (25) 


12.195 (40) 


6 


4 


10^ 


0.26958 (49) 




199.62 (82) 


24.970 (13) 


8 


4 


10« 


0.20015 (66) 




352.6 (1.9) 


41.379 (31) 


10 


4 


10^ 


0.15891 (64) 




606.6 (4.4) 


61.220 (70) 


12 


4 


10^ 


0.13291 (59) 




828.4 (6.9) 


84.45 (10) 


14 


4 


10^ 


0.11277 (61) 




1083 (10) 


110.59 (10) 


16 


4 


10^ 


0.09733 (68) 




1446 (16) 


139.40 (28) 


18 


4 


10^ 


0.08695 (64) 




1846 (24) 


171.92 (27) 


20 


4 


10^ 


0.07852 (58) 




2424 (32) 


206.50 (38) 


30 


4 


10^ 


0.05269 (63) 




5253 (108) 


422.0 (1.2) 


40 


4 


10*^ 


0.03922 (66) 


1.33 (32) 


10325 (281) 


688.6 (2.4) 


50 


4 


10^ 


0.03198 (68) 




16726 (630) 


1021.2 (3.8) 


60 


4 


10^ 


0.02620 (67) 




25525 (1210) 


1412.4 (8.4) 


74 


4 


10^ 


0.02098 (66) 




39769 (2669) 


2016.0 (13.8) 


80 


4 


10^ 


0.02050 (73) 




45501 (3157) 


2355.8 (18.4) 


100 


4 


10^ 


0.01525 (68) 


1.33 (36) 


72877 (6720) 


3389.5 (53.1) 


CXO 


F^ 




0.00033 (34) 








exact 


F^ 













oo 


F^ 




— 


1.33 (26) 






exact 


Fo 




— 


1.764... 
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Table 2: The tensions F£ and effective amplitudes F£ in D 
it = 0.06237). 



2 for /? = 0.47 



L 


No of sweeps 


F[ 


F[ 


T-L 


XL 


10 


4-106 


0.25705 (50) 




629 (5) 


75.123 (47) 


20 


4-106 


0.17959 (68) 


1.94 (67) 


3516 (59) 


292.89 (23) 


40 


4-10^ 


0.14613 (64) 


1.83 (29) 


24267 (1107) 


1164.6 (4) 


60 


4-10^ 


0.13480 (76) 


1.76 (54) 


81338 (6701) 


2614.1 (1.4) 


80 


8-10^ 


0.12969 (60) 


1.96 (17) 


187780 (20224) 


4651.5 (1.4) 


100 


8-10^ 


0.12814 (89) 


1.89 (36) 


571278 (80200) 


7264.9 (2.1) 


oo 


F%m° : 


0.11526 (80) 






0.8522 (1) 


exact 


F^mO : 


0.11492 






0.85196 


oo 


F^-- 


1.848 (12) 


1.92 (17) 






exact 


F^-- 


1.843 


1.93 
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Table 3: The tensions -F£ and effective amplitudes -F£ in D 
(t = 0.11863). 



2 for /5 = 0.5 



L 


No of sweeps 


n 


n 


rL 


XL 


2 


4 


10*^ 


0.6516 (79) 




16.57 (5) 


3.6692 (96) 


4 


4 


106 


0.5852 (54) 




94.81 (36) 


13.907 (61) 


6 


4 


10^ 


0.4780 (06) 




256 (12) 


30.597 (10) 


8 


4 


10^ 


0.4028 (04) 




498 (3) 


53.864 (148) 


10 


4 


10^ 


0.3641 (04) 




856 (43) 


83.792 (32) 


12 


4 


10^ 


0.3369 (14) 




1365 (14) 


120.329 (98) 


14 


4 


10^ 


0.3201 (07) 




1909 (24) 


163.469 (49) 


16 


4 


10^ 


0.3075 (06) 




2851 (45) 


213.330 (72) 


18 


4 


10^ 


0.2991 (04) 




3669 (68) 


269.830 (62) 


20 


4 


10^ 


0.2914 (07) 




5369 (118) 


332.992 (84) 


30 


4 


10^ 


0.2689 (07) 




16688 (550) 


748.44 (18) 


40 


4 


10^ 


0.2580 (05) 




37875 (1861) 


1329.42 (27) 


50 


4 


10*^ 


0.2509 (09) 




80526 (7562) 


2076.85 (60) 


60 


8 


10^ 


0.2478 (07) 




135908 (10444) 


2990.66 (34) 


74 


8 


10^ 


0.2443 (09) 




376843 (32804) 


4548.04 (66) 


80 


8 


10^ 


0.2442 (07) 




585925 (83415) 


5315.71 (70) 


100 


8 


10^ 


0.2403 (12) 


2.04 (7) 


871335 (216745) 


8305.01 (90) 


CXD 


F^mO : 


0.2281(05) 






0.9113 (3) 


exact 


F^m° : 


0.22806... 






0.9113 


oo 


^0 ■ 


1.923 (4) 


2.04 (7) 






exact 


rps . 
^0 ■ 


1.923 


2.10 
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Table 4: The tensions F£ and effective amplitudes F[ in D = 3 as function 
of L in D = 3 for /3 = 0.227 (t = 0.0236). 



L 


No of sweeps 


n 


n 


tl 


XL 


14 


4 


10^ 


0.01146 (10) 




2963 (50) 


639.41 (69) 


16 


4 


10^ 


0.01099 (07) 




4238 (83) 


951.1 (1.3) 


20 


4 


10^ 


0.01092 (08) 


1.291 (20) 


8501 (370) 


1853.8 (2.1) 


24 


4 


10^ 


0.01113 (07) 


1.405 (69) 


14455 (576) 


3188.7 (3.0) 


26 


4 


10^ 


0.01130 (05) 


1.461 (57) 


19626 (907) 


4046.8 (3.4) 


28 


4 


10^ 


0.01149 (04) 


1.465 (71) 


24258 (1142) 


5051 (4) 


32 


8 


10^ 


0.01185 (06) 


1.538 (88) 


36323 (2570) 


7530 (5) 


CXD 


F',m^ : 


0.01293 (17) 






0.4787 (3) 


oo 


i^o: 


1.455 (19) 


1.591 (62) 
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Table 5: The tensions -F£ and effective amplitudes F£ in D 
it = 0.0446). 



3 for /3 = 0.232 



L 


No of sweeps 


n 


n 


tl 


XL 


2 


4-106 


0.16391 (29) 




38.672 (71) 


4.8933 (18) 


4 


4-106 


0.05297 (27) 




160.81 (58) 


26.861 (28) 


6 


4 -106 


0.03403 (15) 




422.1 (2.5) 


78.00 (8) 


8 


4-106 


0.02779 (13) 




880.7 (7.6) 


176.15 (18) 


10 


4-106 


0.02561 (11) 




1566 (19) 


340.76 (29) 


12 


4-106 


0.02485 (12) 




2671 (41) 


587.10 (41) 


14 


4-106 


0.02479 (10) 




4168 (82) 


929.26 (50) 


16 


4-106 


0.02521 (11) 




6604 (151) 


1383.2 (9) 


20 


4.2 - 106 


0.02630 (11) 




13269 (546) 


2698.1 (1.8) 


24 


8-106 


0.02752 (05) 


1.579 (20) 


23528 (1068) 


4652.4 (1.4) 


26 


4-106 


0.02794 (07) 


1.613 (41) 


35032 (1985) 


5912.8 (2.7) 


28 


4-106 


0.02847 (07) 


1.591 (39) 


49670 (4342) 


7383.7 (3.4) 


32 


6-106 


0.02906 (08) 


1.581 (75) 


91496 (6902) 


11011.8 (5.3) 


oo 


F^m° : 


0.03140 (14) 






0.5797 (1) 


oo 


F' : 


1.581 (07) 


1.616 (82) 
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Table 6: The tensions Fl and effective amplitudes Fl in D 
it = 0.0912). 



3 for /5 = 0.2439 



L 


No of sweeps 


n 


n 


tl 


XL 


2 


4-106 


0.19535 (30) 




38.412 (51) 


5.2022 (22) 


4 


4-106 


0.08449 (21) 




193.8 (8) 


33.592 (23) 


6 


4 -106 


0.06691 (19) 




563 (4) 


110.058 (77) 


8 


4-106 


0.06327 (16) 




1221 (12) 


259.85 (10) 


10 


4-106 


0.06397 (13) 




2396 (35) 


505.738 (192) 


12 


4-106 


0.06556 (14) 




4339 (81) 


871.35 (40) 


14 


4-106 


0.06785 (12) 




7433 (196) 


1382.61 (50) 


16 


3.6 - 106 


0.06989 (13) 




12970 (494) 


2061.3 (6) 


20 


4-106 


0.07254 (09) 


1.694 (31) 


40374 (2296) 


4025.7 (1.5) 


24 


8-106 


0.07349 (11) 


1.638 (53) 


76965 (7056) 


6954 (1) 


28 


4-106 


0.07363 (14) 


1.599 (44) 


238560 (29158) 


11042 (3) 


32 


6.6 - 106 


0.07367 (10) 


1.687 (33) 


508857 (94588) 


16480 (4) 


oo 


F^mO : 


0.07403 (30) 






0.70913 (7) 


oo 


-^0 • 


1.513 (06) 


1.647 (59) 
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Table 7: Cusp estimates m^. 



D = 2 


/3 = 0.5 


D = 3 


/3 = 0.2439 


L 


ml 


L 


ml 


60 


0.346 (11) 


24 


0.327 (07) 


80 


0.351 (21) 


28 


0.336 (14) 


100 


0.351 (45) 


32 


0.332 (14) 



Figure Captions 

Figure 1: Probability distribution for the magnetization on the L = 32 
lattice for (3 = 0.227 in the 3D Ising model. This distribution was obtained 
using the multimagnetic parameters as estimated by our finite size scaling 
prediction. 

Figure 2: 2D Magnetic probability density Pl{M) for P = 0.5. 

Figure 3: 2D effective tensions as a function of the lattice size for f3 = f3c, 
0.47, 0.5. The curves are our best fits as described in the text. 

Figure 4: 2D canonical heatbath and multimagnetical tunneling times 
versus lattice size at /? = /3c. The lines are power law fits and the solid line 
correspond to the multimagnetical simulation. 

Figure 5: 3D Magnetic probability density PLiM) for p = 0.2439. 

Figure 6: 3-D effective tensions as a function of the lattice size for 
P = 0.227, 0.232, 0.2439. The curves are our asymptotic fits as described in 

the text. 

Figure 7: Typical m ^ configuration at /? = 0.5 on a 100^ lattice. 

Figure 8: Relying on 2D Ising model analytical results 0: Equilibrium 
shapes of droplets for several temperatures or /3- values. r{9) denotes the 
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angle dependent radius of the droplet. 

Figure 9: Numerical estimate of the cusp m'^ at /? = 0.5 in 2D. The 
dashed line indicates a possible extrapolation to the theoretically known in- 
finite volume value (solid circle). 

Figure 10: We plot 0.5 ■ ilj{m) at (3 = 0.5 in 2D. The curve corresponds 
to a fit according to eq. (30). 

Figure 11: Surface area per link versus magnetization at /9 = 0.2439 in 
3D. 
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